Quantum Monte Carlo simulation of a two-dimensional Bose gas 
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The equation of state of a homogeneous two-dimensional Bose gas is calculated using quantum 
Monte Carlo methods. The low-density universal behavior is investigated using different interatomic 
model potentials, both finite- ranged and strictly repulsive and zero-ranged supporting a bound state. 
The condensate fraction and the pair distribution function are calculated as a function of the gas 
parameter, ranging from the dilute to the strongly correlated regime. In the case of the zero-range 
pseudopotential we discuss the stability of the gas-like state for large values of the two-dimensional 
scattering length, and we calculate the critical density where the system becomes unstable against 
cluster formation. 

PACS numbers: 



I. INTRODUCTION 

Degenerate low-dimensional gases are presently at- 
tracting considerable interest as model systems to inves- 
tigate beyond mean-field effects and phenomena where 
many-body correlations, thermal and quantum fluctua- 
tions play a relevant role In two dimensions (2D), it 
is well known that thermal excitations destroy long-range 
order and Bose-Einstein condensation (BEC) can not ex- 
ist in large bosonic systems at finite temperature 0- 
Nevertheless, a defect-mediated phase transition from 
a high-temperature normal fluid to a low-temperature 
superfluid was predicted by Berezinskii, Kosterlitz and 
Thouless (BKT) and has been observed in thin films 
of hquid '^He ^. At zero temperature, fluctuations arc 
instead negligible and BEC is possible. The ground-state 
energy per particle of a homogeneous dilute Bose gas in 
2D has been first calculated by Schick ^ and is given by 
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where m is the mass of the particles, n is the number den- 
sity and a2D is the 2D scattering length. The above result 



holds for small values of the gas parameter, 



^2D 



< 1, 



and can be obtained within a mean-field approximation 
using the coupling constant 
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The presence at low energy of a density dependent cou- 
pling constant is a peculiar feature of 2D systems. 

Bose gases in quasi-2D have been realized with spin- 
polarized hydrogen on a liquid-helium surface and 
with alkali-metal atoms confined in highly anisotropic 
harmonic traps 0; B I3- these systems, the trans- 
verse motion of the atoms is frozen to zero-point oscil- 
lations resulting in fully 2D kinematics, but the inter- 
atomic scattering processes are still three-dimensional 



(3D). In fact, the 3D s-wave scattering length a^u is 
always much smaller than the characteristic length of 
the tight transverse confinement. Theoretical studies 
of quasi-2D trapped systems predict the existence of a 
low-temperature quasicondensate phase [Tol [Tll | whose 
properties have been investigated using mean-field ap- 
proaches The nature of the transition in confined 
systems and whether it belongs to the BEC or BKT uni- 
versality class is still an open problem both experimen- 
tally and theoretically. 

An important theoretical prediction for quasi-2D Bose 
gases in harmonic traps is the renormalization of the 
coupling constant due to the tight confinement in one 
direction. To logarithmic accuracy the result is given 
by [Hill 
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where = \JfilmiOz is the harmonic oscillator length, 
fixed by the frequency lo^ of the tight transverse confine- 
ment, and n(0) is the 2D density in the center of the trap. 
If a3£) ^ , one can neglect the log term in Eq. and 
the resulting coupling constant is determined by the 3D 
scattering length. On the contrary, if a-^D 3> a^, g2o 
becomes independent of the value of aj^u and a regime 
of pure 2D scattering is achieved with a20 = Oz. This 
universal regime, where the properties of the system only 
depend on density and not on interaction, is highly in- 
teresting and can be realized in trapped gases by using a 
Feshbach resonance j3 to achieve large values for a-^u- 
In this work we investigate the ground-state proper- 
ties of a strictly 2D homogeneous Bose gas using quan- 
tum Monte Carlo techniques. We determine the equation 
of state over a wide range of values of the gas parame- 
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ter from the extremely dilute regime 
the strongly correlated regime, na^j^ ~ 0.1. The calcu- 
lations of the energy per particle are carried out using 
three different interatomic model potentials: two repul- 
sive finite-ranged potentials and a zero-range pseudopo- 
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tential which supports a two-body bound state. We in- 
vestigate beyond mean-field corrections to the equation 
of state and their dependence on the gas parame- 
ter na2]2)- In the case of the zero-range pseudopotential, 
from a study of the system compressibihty as a function 
of the gas parameter we estimate that na^jj — 0.04 is 
the critical density at which the gas-like state becomes 
unstable against cluster formation. This result puts an 
upper limit to the value of density that can be reached 
in harmonic traps when one enters the regime of pure 2D 
scattering. We also report results of the pair distribution 
function and the condensate fraction as a function of the 
gas parameter. 

The structure of the paper is as follows. In Section UTI 
we introduce the many-body Hamiltonian and the model 
interatomic potentials we use in the simulations. We also 
discuss the quantum Monte Carlo methods employed in 
the present study and the choice of the trial wave func- 
tions. Section Hill containcs the results obtained and Sec- 
tion Uvl draws our conclusions. 



II. METHOD 

We consider a homogeneous system of N spinless 
bosons in 2D described by the many-body Hamiltonian 
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where r.^ — Xii + denotes the 2D position vector of 
the i-th particle and V{r) is the two-body interatomic 
potential. We use different choices for the interatomic 
potential V{r): 

(i) Hard-disk (HD) potential defined by 



oo (r < a2D) 
(r > a2D) , 



(5) 



in which case the hard-disk diameter corresponds to the 
scattering length a2D- 

(ii) Soft-disk (SD) potential defined by 



Vo {r<R) 
(r > i?) 



(6) 



in which case the scattering length is given in terms of the 
modified Bessel function Io{x) and its derivative Io{x), 



0-2D = i?exp 



1 loiRKo) 



RK^I'^iRKo) 



(7) 



where — mVo/h^ with Vb > 0. For a finite Vb one has 
always R > a2D\ if Vb +oo the SD and HD potentials 
coincide with a2D — R- In the present calculations, the 
range of the SD potential is kept fixed to the value R — 
5a2D and the height Vq is determined through Eq. {Tj) to 
give the desired value of a2D- 



(iii) Zero-range pseudopotential (PP) defined by ^i 



m log(ga2D) 
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(8) 



where the second term in the square brackets is a regu- 
larizing operator with q an arbitrary wave vector. The 
potential V^^ (r) supports a two-body bound state with 
energy eb — —^h?/{'ma'^jje'^'') and wave function /fc(r) = 
Ko{2r / {a2DG')), where Ki^{x) is the modified Bessel 
function and 7 ~ 0.577 is the Euler constant. 

The calculations with the repulsive potentials V^^ (r) 
and V^^{r) are carried out using the diffusion Monte 
Carlo (DMC) method. This technique solves the many- 
body time-independent Schrodinger equation by evolving 
the function /(R, t) ~ V-'t(R-)^(R-, t) in imaginary time 
T = it/h according to the time-dependent Schrodinger 
equation 



g/(R,r) 
dr 
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Here 'I'(R, t) denotes the wave function of the sys- 
tem and f/jTlR) is a trial function used for importance 
sampling. In the above equation R = (ri, r^r), 
Sl(R) = V'T(R)"^^fV'T(R) denotes the local energy, 
F(R) = 27/it(R)"^Vri/'t(R) is the quantum drift force, 
D — h^/{2m) plays the role of an effective diffusion 
constant, and Eref is a reference energy introduced to 
stabilize the numerics. The energy and other observ- 
ables of the ground state of the system are calculated 
from averages over the asymptotic distribution function 
/(R, r 00). Apart from statistical errors, the DMC 
method determines the ground-state energy of a system 
of N bosons exactly. 

In the case of the pseudopotential V^^{r), Eq. 0, 
the gas-like state is not the ground state of the system, 
but is a highly-excited state of the Hamiltonian. We ex- 
pect, however, that for small values of the gas parameter, 
na^jj ^ 1, the gas-like state is stable against formation 
of many-body bound states. The calculations with the 
PP potential are performed at the level of the variational 
Monte Carlo (VMC) method to avoid the technical dif- 
ficulties that arise in the study of excited states in the 
DMC method. Nevertheless, as it will be shown in Sec- 
tion the accuracy achieved by using this variational 
approach is very high. The VMC technique is based on 
the variational principle 



(V>T|g|^T) 
{lpT\i>T) 



> E 



(10) 



which is here applied to excited states of the Hamiltonian 
H. For a trial wave function ipT with a given symmetry, 
the variational estimate provides an upper bound to the 
energy E of the lowest excited state of the Hamiltonian 
H with that symmetry. 
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In all the calculations the trial wave function '0t(R-) is 
taken of the Jastrow form 

V'T(ri,...,rA,)=n/(r,,), (11) 

where the correlation factor f{r) is chosen as the solu- 
tion of the two-body Schrodinger equation subject to the 
boundary conditions f{r > d) — 1 and f'{r > d) — 
on the wave function and its first derivative respectively, 
with d a variational parameter. For the repulsive poten- 
tials, HD and SD, fir) corresponds to the ground-state 
two-body wave function. On the contrary, in the case of 
the PP interaction, /(r) is chosen as the two-body wave 
function corresponding to the lowest excited state with 
positive energy. The pseudopotential V^^ imposes the 
following boundary condition on the pair wave function 
for zero interparticle distance 

[f{r)-\og{qr)rf'{T%^^ \og{qa2D) ' ^ ' 

The free two-particle problem in 2D — ?i^/mV^/(r) = 
Tt^k^ /mf{r), holding for r 7^ 0, gives the general solution 

/(r) = AJo(fcr)+Bro(fcr) , (13) 

in terms of the Jo (a;) and lo(x) Bessel functions. The 
boundary condition (|12|) is imposed by the relation 

a,^ = ^exp(^-— j . (14) 

In addition we require that /(r) possess only one node for 
distances < r < d. The boundary conditions at r = d 
together with (|14|l fix the values of the constants A and 
B and of the wavevector k. If d ^ 02,0 , the eigenenergy 
f?k'^ jm is small and the node of /(r) is located at r = 

0.20- 

All the Monte Carlo calculations have been carried out 
considering periodic boundary conditions with a cubic 
simulation box whose size L is fixed by the density n and 
by the total number of atoms: n = N/L"^. The healing 
length d is considered everywhere d = L/2 since a VMC 
optimization has shown that this is a good option. 

III. RESULTS 

The numerical simulations are carried out with N = 
512 particles for the lowest densities {na^jj < lO""*) and 
with N — 256 particles for the rest. The DMC results 
for the energy per particle corresponding to the HD po- 
tential are reported in Table HI as a function of the gas 
parameter. Variational results on the same system have 
been obtained in the dilute regime by Mazzanti et al. [l^ 
using Correlated Basis Functions theory and at high den- 
sities by Lei Xing {V^ and Soli's [13 using VMC. In Fig.Hl 
the DMC results for the HD system are compared with 



TABLE I: Energy per particle E/N corresponding to the 
model potentials V^'^ir) (DMC results), V'^^ir) (DMC re- 
sults) and y^^(r) (VMC results) [energies are in units of 
?iV(2maiB)]. 
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FIG. 1: Energy per particle as a function of the gas parameter 
for the HD potential. Solid line: mean-field result, Eq. Q. 
Error bars are smaller than the size of the symbols. 



the mean- field result, Eq. We find a remarkably 

good agreement up to large values of na\jj. It is worth 
noticing that the range of validity of the mean-field ap- 
proximation is in 2D significantly larger than the one 
observed in 3D Beyond mean- field corrections to 

the equation of state are shown in Fig. |21 where we plot 
the difference {E — Emf)/N for the three model poten- 
tials: V"^{r) (DMC results), V^^{t) (DMC results) 
and y^^(r) (VMC results). Focusing on the exact DMC 
results, one can notice the universal behavior of the equa- 
tion of state for na^jj <C 1. The results of the HD and SD 
model potentials are practically indistinguishable up to 
na2u ~ 10~^ (see Table QJ, showing that beyond mean- 
field corrections to Eq. only depend on the value 
of the gas parameter na^jj and not on the details of the 
interatomic potential. By increasing na^D the results ob- 
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FIG. 2: Beyond mean-field corrections to tlie equation of 
state. Circles, HD potential; triangles: SD potential; squares: 
PP potential. Dashed line: numerical fit to the results of 
the HD potential, Eq. 11611 . Error bars are smaller than the 
size of the symbols. Inset: region of extremely dilute systems 




FIG. 3: Energy per particle and inverse compressibility as 
a function of the gas parameter. Symbols are as in Fig. |21 
Thick and thin dashed lines: best fit p6|l to the HD equation 
of state and corresponding mc^, respectively. Thick and thin 
solid line: polynomial best fit to the PP equation of state and 
corresponding mc^, respectively. Error bars are smaller than 
the size of the symbols. 



tained with the SD potential (triangles in Fig. 12) start to 
deviate from the ones corresponding to the HD potential 
(circles in Fig.[21l approximately about na^jj ~ 0.01. The 
VMC results for the PP potential are upper bounds to 
the eigenenergy of the gas-like state. However, looking at 
TablelHand Fig.|21 one can see that they are nearly indis- 
tinguishable from the exact HD and SD energies in the 
low-density regime. A similar accuracy for the energy of 
the gas-like state is expected also for larger values of the 
gas parameter, where the results of the PP interaction 
start to separate from the HD ones as the effects of the 
finite interaction range become more important. It is also 
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FIG. 4: Condensate fraction as a function of the gas parame- 
ter. Circles, HD potential; triangles, SD potential. Solid line: 
result from Bogoliubov theory, Eq. II17II . 
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FIG. 5: Pair distribution function of a HD gas obtained for 



three values of the gas parameter: na2£) 



10" 



solid line; 



: 0.01, dashed line; noon = 0.1, dotted line. 



worth noticing that for the smallest values of na^jj (inset 
of Fig. 121) the beyond mean-field correction is negative. 
This is consistent with the perturbation expansion po[ 



E-Emf _ 27r?i^n log[log(l/na^£,)] 



^2D 
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(15) 



TV TO [log(r7,a2 

which is expected to hold if log[log(l/na2£))] ^ 1. For 
our smallest value of the gas parameter, na^jj — 10~^, 
one finds log[log(l/na2£))] ~ 2.8. We fit the HD equation 
of state using the expression 



E — Emf 
N 



(16) 



-A 



log[log(l/na2^) 



B 



[\og{r 



''2D 



[log(r 



'■2D 



where A and B are free parameters. The best fit to all 
HD points, with the exclusion of the one at the highest 
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density na^^, = 0.1, gives A = 0.86 and B = 2.26 with 
X^/v — 1.2. The resulting curve is shown in Fig.|21with 
a dashed hne. We notice that the value A = 0.86 is close 
to the prediction A — 1 from the expansion (|15|l . 

In Fig. we show an enlargement of the equation of 
state in the gas parameter range 0.01 < na^u < 0.1. 
In this region, the VMC results corresponding to the 
V^^{r) potential significantly deviate from the HD equa- 
tion of state (thick dashed line). By using a polyno- 
mial fit to the PP equation of state (thick solid line), 
we calculate the inverse compressibility of the system 
mc^ = nd^/dn, where c is the speed of sound and 
/i — dE/dN is the chemical potential. The inverse com- 
pressibility (thin solid line) exhibits a maximum and then 
drops to zero for na\jj ~ 0.04. For comparison, the in- 
verse compressibility of the HD gas as obtained from the 
fit (|16|l is shown with a thin dashed line. We interpret 
the vanishing of mc^ as the onset of instability against 
cluster formation. 

The results for the condensate fraction corresponding 
to the HD and SD potentials are shown in Fig. 01 as 
a function of the gas parameter. The condensate frac- 
tion, Nq/N, is obtained from the long-range behavior of 
the one-body density matrix Nq/N = Huir^oo pir) (see 
Ref. for further details). The Bogoliubov theory ap- 
plied to the 2D Bose gas gives the result Q 

^.1 + I (17) 

which should be contrasted with the 3D result Nq/N = 
l-8Vna|^/(30F). From Fig. H one sees that for small 
values of the gas parameter the results of both the HD 
and SD potential agree with Eq. H17|) . For larger val- 
ues of na|^ the Bogoliubov result H17() overestimates the 
condensate fraction of the HD gas. In the case of the SD 
potential we notice that A'o /N first decreases by increas- 
ing the gas parameter and then increases. This happens 



at densities where the interaction ranges of the particles 
start to overlap resulting in a reduction of interaction 
effects. It is worth noticing that in the region where Bo- 
goliubov theory applies, for the same value of the gas 
parameter the condensate depletion is larger in 2D than 
in 3D, showing that quantum fluctuations are more ef- 
fective in systems with reduced dimensionality. 

Finally in Fig. El we report results of the pair distribu- 
tion function g2 (r) corresponding to the HD potential for 
three different values of the gas parameter na^jj = 10~^, 
0.01, 0.1. The function (72(f) gives the probability that 
two particles will be separated by a distance r. A shell- 
like structure in the pair distribution function is visible 
only for the highest density. The reduction from three 
to two dimensions changes the long-range behavior of 
g2{r) turning the dependence {g2^{r) — 1) oc into 
{g2^{r) — 1) (X and therefore the relevance of long- 
range correlations is enhanced in 2D. 

IV. CONCLUSIONS 

We have carried out a study of the ground-state prop- 
erties of a 2D homogeneous Bose gas with quantum 
Monte Carlo techniques. The universal behavior of the 
equation of state for small values of the gas parameter 
has been investigated by analyzing both mean-field and 
beyond mean-field contributions. By using a zero-range 
pseudopotential to describe interatomic interactions we 
have estimated that na^jj — 0.04 is the critical density 
at which the system becomes unstable against cluster 
formation. These results are relevant for present and 
future experiments with ultracold Bose gases in highly 
anisotropic harmonic traps. 
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